home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 November / macformat-018.iso / Utility Spectacular / Utilities / Calc / zmath.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-10-10  |  31.5 KB  |  1,777 lines  |  [TEXT/????]
  1. /*
  2.  * Copyright (c) 1992 David I. Bell
  3.  * Permission is granted to use, distribute, or modify this source,
  4.  * provided that this copyright notice remains intact.
  5.  *
  6.  * Extended precision integral arithmetic primitives
  7.  */
  8.  
  9. #include <stdio.h>
  10. #include "xmath.h"
  11.  
  12.  
  13. HALF _twoval_[] = { 2 };
  14. HALF _zeroval_[] = { 0 };
  15. HALF _oneval_[] = { 1 };
  16. HALF _tenval_[] = { 10 };
  17.  
  18. ZVALUE _zero_ = { _zeroval_, 1, 0};
  19. ZVALUE _one_ = { _oneval_, 1, 0 };
  20. ZVALUE _ten_ = { _tenval_, 1, 0 };
  21.  
  22. /*
  23.  * mask of given bits, rotated thru all bit positions twice
  24.  *
  25.  * bitmask[i]      (1 << (i-1)),  for  -BASEB*4<=i<=BASEB*4
  26.  * rotmask[j][k] (1 << ((j+k-1)%BASEB)),  for  -BASEB*2<=j,k<=BASEB*2
  27.  */
  28. static HALF *bmask;        /* actual rotation thru 8 cycles */
  29. static HALF **rmask;        /* actual rotation pointers thru 2 cycles */
  30. HALF *bitmask;            /* bit rotation, norm 0 */
  31. #if 0
  32. static HALF **rotmask;        /* pointers to bit rotations, norm index */
  33. #endif
  34.  
  35. BOOL _math_abort_;        /* nonzero to abort calculations */
  36.  
  37.  
  38. static char *abortmsg = "Calculation aborted";
  39. static char *memmsg = "Not enough memory";
  40.  
  41. static void dadd proto((ZVALUE z1, ZVALUE z2, long y, long n));
  42. static BOOL dsub proto((ZVALUE z1, ZVALUE z2, long y, long n));
  43. static void dmul proto((ZVALUE z, FULL x, ZVALUE *dest));
  44.  
  45.  
  46. #ifdef ALLOCTEST
  47. static long nalloc = 0;
  48. static long nfree = 0;
  49. #endif
  50.  
  51.  
  52. HALF *
  53. alloc(len)
  54.     long len;
  55. {
  56.     HALF *hp;
  57.  
  58.     if (_math_abort_)
  59.         error(abortmsg);
  60.     hp = (HALF *) malloc(len * sizeof(HALF) + 2);
  61.     if (hp == 0)
  62.         error(memmsg);
  63. #ifdef ALLOCTEST
  64.     ++nalloc;
  65. #endif
  66.     return hp;
  67. }
  68.  
  69. #ifdef ALLOCTEST
  70. void
  71. freeh(h)
  72.     HALF *h;
  73. {
  74.     if ((h != _zeroval_) && (h != _oneval_)) {
  75.         free(h);
  76.         ++nfree;
  77.     }
  78. }
  79.  
  80.  
  81. void
  82. allocStat()
  83. {
  84.     fprintf(stderr, "nalloc: %ld nfree: %ld kept: %ld\n",
  85.         nalloc, nfree, nalloc - nfree);
  86. }
  87. #endif
  88.  
  89.  
  90. /*
  91.  * Convert a normal integer to a number.
  92.  */
  93. void
  94. itoz(i, res)
  95.     long i;
  96.     ZVALUE *res;
  97. {
  98.     long diddle, len;
  99.  
  100.     res->len = 1;
  101.     res->sign = 0;
  102.     diddle = 0;
  103.     if (i == 0) {
  104.         res->v = _zeroval_;
  105.         return;
  106.     }
  107.     if (i < 0) {
  108.         res->sign = 1;
  109.         i = -i;
  110.         if (i < 0) {    /* fix most negative number */
  111.             diddle = 1;
  112.             i--;
  113.         }
  114.     }
  115.     if (i == 1) {
  116.         res->v = _oneval_;
  117.         return;
  118.     }
  119.     len = 1 + (((FULL) i) >= BASE);
  120.     res->len = len;
  121.     res->v = alloc(len);
  122.     res->v[0] = (HALF) (i + diddle);
  123.     if (len == 2)
  124.         res->v[1] = (HALF) (i / BASE);
  125. }
  126.  
  127.  
  128. /*
  129.  * Convert a number to a normal integer, as far as possible.
  130.  * If the number is out of range, the largest number is returned.
  131.  */
  132. long
  133. ztoi(z)
  134.     ZVALUE z;
  135. {
  136.     long i;
  137.  
  138.     if (isbig(z)) {
  139.         i = MAXFULL;
  140.         return (z.sign ? -i : i);
  141.     }
  142.     i = (istiny(z) ? z1tol(z) : z2tol(z));
  143.     return (z.sign ? -i : i);
  144. }
  145.  
  146.  
  147. /*
  148.  * Make a copy of an integer value
  149.  */
  150. void
  151. zcopy(z, res)
  152.     ZVALUE z, *res;
  153. {
  154.     res->sign = z.sign;
  155.     res->len = z.len;
  156.     if (isleone(z)) {    /* zero or plus or minus one are easy */
  157.         res->v = (z.v[0] ? _oneval_ : _zeroval_);
  158.         return;
  159.     }
  160.     res->v = alloc(z.len);
  161.     copyval(z, *res);
  162. }
  163.  
  164.  
  165. /*
  166.  * Add together two integers.
  167.  */
  168. void
  169. zadd(z1, z2, res)
  170.     ZVALUE z1, z2, *res;
  171. {
  172.     ZVALUE dest;
  173.     HALF *p1, *p2, *pd;
  174.     long len;
  175.     FULL carry;
  176.     SIUNION sival;
  177.  
  178.     if (z1.sign && !z2.sign) {
  179.         z1.sign = 0;
  180.         zsub(z2, z1, res);
  181.         return;
  182.     }
  183.     if (z2.sign && !z1.sign) {
  184.         z2.sign = 0;
  185.         zsub(z1, z2, res);
  186.         return;
  187.     }
  188.     if (z2.len > z1.len) {
  189.         pd = z1.v; z1.v = z2.v; z2.v = pd;
  190.         len = z1.len; z1.len = z2.len; z2.len = len;
  191.     }
  192.     dest.len = z1.len + 1;
  193.     dest.v = alloc(dest.len);
  194.     dest.sign = z1.sign;
  195.     carry = 0;
  196.     pd = dest.v;
  197.     p1 = z1.v;
  198.     p2 = z2.v;
  199.     len = z2.len;
  200.     while (len--) {
  201.         sival.ivalue = ((FULL) *p1++) + ((FULL) *p2++) + carry;
  202.         *pd++ = sival.silow;
  203.         carry = sival.sihigh;
  204.     }
  205.     len = z1.len - z2.len;
  206.     while (len--) {
  207.         sival.ivalue = ((FULL) *p1++) + carry;
  208.         *pd++ = sival.silow;
  209.         carry = sival.sihigh;
  210.     }
  211.     *pd = (HALF)carry;
  212.     quicktrim(dest);
  213.     *res = dest;
  214. }
  215.  
  216.  
  217. /*
  218.  * Subtract two integers.
  219.  */
  220. void
  221. zsub(z1, z2, res)
  222.     ZVALUE z1, z2, *res;
  223. {
  224.     register HALF *h1, *h2, *hd;
  225.     long len1, len2;
  226.     FULL carry;
  227.     SIUNION sival;
  228.     ZVALUE dest;
  229.  
  230.     if (z1.sign != z2.sign) {
  231.         z2.sign = z1.sign;
  232.         zadd(z1, z2, res);
  233.         return;
  234.     }
  235.     len1 = z1.len;
  236.     len2 = z2.len;
  237.     if (len1 == len2) {
  238.         h1 = z1.v + len1 - 1;
  239.         h2 = z2.v + len2 - 1;
  240.         while ((len1 > 0) && ((FULL)*h1 == (FULL)*h2)) {
  241.             len1--;
  242.             h1--;
  243.             h2--;
  244.         }
  245.         if (len1 == 0) {
  246.             *res = _zero_;
  247.             return;
  248.         }
  249.         len2 = len1;
  250.     }
  251.     dest.sign = z1.sign;
  252.     carry = ((len1 < len2) || ((len1 == len2) && ((FULL)*h1 < (FULL)*h2)));
  253.     h1 = z1.v;
  254.     h2 = z2.v;
  255.     if (carry) {
  256.         carry = len1;
  257.         len1 = len2;
  258.         len2 = carry;
  259.         h1 = z2.v;
  260.         h2 = z1.v;
  261.         dest.sign = !dest.sign;
  262.     }
  263.     hd = alloc(len1);
  264.     dest.v = hd;
  265.     dest.len = len1;
  266.     len1 -= len2;
  267.     carry = 0;
  268.     while (--len2 >= 0) {
  269.         sival.ivalue = (BASE1 - ((FULL) *h1++)) + *h2++ + carry;
  270.         *hd++ = BASE1 - sival.silow;
  271.         carry = sival.sihigh;
  272.     }
  273.     while (--len1 >= 0) {
  274.         sival.ivalue = (BASE1 - ((FULL) *h1++)) + carry;
  275.         *hd++ = BASE1 - sival.silow;
  276.         carry = sival.sihigh;
  277.     }
  278.     if (hd[-1] == 0)
  279.         trim(&dest);
  280.     *res = dest;
  281. }
  282.  
  283.  
  284. #if 0
  285. /*
  286.  * Multiply two integers together.
  287.  */
  288. void
  289. zmul(z1, z2, res)
  290.     ZVALUE z1, z2, *res;
  291. {
  292.     register HALF *s1, *s2, *sd, *h1;
  293.     FULL mul, carry;
  294.     long len1, len2;
  295.     SIUNION sival;
  296.     ZVALUE dest;
  297.  
  298.     if (iszero(z1) || iszero(z2)) {
  299.         *res = _zero_;
  300.         return;
  301.     }
  302.     if (isone(z1)) {
  303.         zcopy(z2, res);
  304.         return;
  305.     }
  306.     if (isone(z2)) {
  307.         zcopy(z1, res);
  308.         return;
  309.     }
  310.     dest.len = z1.len + z2.len;
  311.     dest.v = alloc(dest.len);
  312.     dest.sign = (z1.sign != z2.sign);
  313.     clearval(dest);
  314.     /*
  315.      * Swap the numbers if necessary to make the second one smaller.
  316.      */
  317.     if (z1.len < z2.len) {
  318.         len1 = z1.len;
  319.         z1.len = z2.len;
  320.         z2.len = len1;
  321.         s1 = z1.v;
  322.         z1.v = z2.v;
  323.         z2.v = s1;
  324.     }
  325.     /*
  326.      * Multiply the first number by each 'digit' of the second number
  327.      * and add the result to the total.
  328.      */
  329.     sd = dest.v;
  330.     s2 = z2.v;
  331.     for (len2 = z2.len; len2--; sd++) {
  332.         mul = (FULL) *s2++;
  333.         if (mul == 0)
  334.             continue;
  335.         h1 = sd;
  336.         s1 = z1.v;
  337.         carry = 0;
  338.         len1 = z1.len;
  339.         while (len1 >= 4) {    /* expand the loop some */
  340.             len1 -= 4;
  341.             sival.ivalue = (mul * ((FULL) *s1++)) + ((FULL) *h1) + carry;
  342.             *h1++ = sival.silow;
  343.             sival.ivalue = (mul * ((FULL) *s1++)) + ((FULL) *h1) + ((FULL) sival.sihigh);
  344.             *h1++ = sival.silow;
  345.             sival.ivalue = (mul * ((FULL) *s1++)) + ((FULL) *h1) + ((FULL) sival.sihigh);
  346.             *h1++ = sival.silow;
  347.             sival.ivalue = (mul * ((FULL) *s1++)) + ((FULL) *h1) + ((FULL) sival.sihigh);
  348.             *h1++ = sival.silow;
  349.             carry = sival.sihigh;
  350.         }
  351.         while (--len1 >= 0) {
  352.             sival.ivalue = (mul * ((FULL) *s1++)) + ((FULL) *h1) + carry;
  353.             *h1++ = sival.silow;
  354.             carry = sival.sihigh;
  355.         }
  356.         while (carry) {
  357.             sival.ivalue = ((FULL) *h1) + carry;
  358.             *h1++ = sival.silow;
  359.             carry = sival.sihigh;
  360.         }
  361.     }
  362.     trim(&dest);
  363.     *res = dest;
  364. }
  365. #endif
  366.  
  367.  
  368. /*
  369.  * Multiply an integer by a small number.
  370.  */
  371. void
  372. zmuli(z, n, res)
  373.     ZVALUE z;
  374.     long n;
  375.     ZVALUE *res;
  376. {
  377.     register HALF *h1, *sd;
  378.     FULL low;
  379.     FULL high;
  380.     FULL carry;
  381.     long len;
  382.     SIUNION sival;
  383.     ZVALUE dest;
  384.  
  385.     if ((n == 0) || iszero(z)) {
  386.         *res = _zero_;
  387.         return;
  388.     }
  389.     if (n < 0) {
  390.         n = -n;
  391.         z.sign = !z.sign;
  392.     }
  393.     if (n == 1) {
  394.         zcopy(z, res);
  395.         return;
  396.     }
  397.     low = ((FULL) n) & BASE1;
  398.     high = ((FULL) n) >> BASEB;
  399.     dest.len = z.len + 2;
  400.     dest.v = alloc(dest.len);
  401.     dest.sign = z.sign;
  402.     /*
  403.      * Multiply by the low digit.
  404.      */
  405.     h1 = z.v;
  406.     sd = dest.v;
  407.     len = z.len;
  408.     carry = 0;
  409.     while (len--) {
  410.         sival.ivalue = ((FULL) *h1++) * low + carry;
  411.         *sd++ = sival.silow;
  412.         carry = sival.sihigh;
  413.     }
  414.     *sd = (HALF)carry;
  415.     /*
  416.      * If there was only one digit, then we are all done except
  417.      * for trimming the number if there was no last carry.
  418.      */
  419.     if (high == 0) {
  420.         dest.len--;
  421.         if (carry == 0)
  422.             dest.len--;
  423.         *res = dest;
  424.         return;
  425.     }
  426.     /*
  427.      * Need to multiply by the high digit and add it into the
  428.      * previous value.  Clear the final word of rubbish first.
  429.      */
  430.     *(++sd) = 0;
  431.     h1 = z.v;
  432.     sd = dest.v + 1;
  433.     len = z.len;
  434.     carry = 0;
  435.     while (len--) {
  436.         sival.ivalue = ((FULL) *h1++) * high + ((FULL) *sd) + carry;
  437.         *sd++ = sival.silow;
  438.         carry = sival.sihigh;
  439.     }
  440.     *sd = (HALF)carry;
  441.     quicktrim(dest);
  442.     *res = dest;
  443. }
  444.  
  445.  
  446. /*
  447.  * Divide two numbers by their greatest common divisor.
  448.  * This is useful for reducing the numerator and denominator of
  449.  * a fraction to its lowest terms.
  450.  */
  451. void
  452. zreduce(z1, z2, z1res, z2res)
  453.     ZVALUE z1, z2, *z1res, *z2res;
  454. {
  455.     ZVALUE tmp;
  456.  
  457.     if (isleone(z1) || isleone(z2))
  458.         tmp = _one_;
  459.     else
  460.         zgcd(z1, z2, &tmp);
  461.     if (isunit(tmp)) {
  462.         zcopy(z1, z1res);
  463.         zcopy(z2, z2res);
  464.     } else {
  465.         zquo(z1, tmp, z1res);
  466.         zquo(z2, tmp, z2res);
  467.     }
  468.     freeh(tmp.v);
  469. }
  470.  
  471.  
  472. /*
  473.  * Divide two numbers to obtain a quotient and remainder.
  474.  * This algorithm is taken from
  475.  * Knuth, The Art of Computer Programming, vol 2: Seminumerical Algorithms.
  476.  * Slight modifications were made to speed this mess up.
  477.  */
  478. void
  479. zdiv(z1, z2, res, rem)
  480.     ZVALUE z1, z2, *res, *rem;
  481. {
  482.     long i, j, k;
  483.     register HALF *q, *pp;
  484.     SIUNION pair;        /* pair of halfword values */
  485.     HALF h2, v2;
  486.     long y;
  487.     FULL x;
  488.     ZVALUE ztmp1, ztmp2, ztmp3, quo;
  489.  
  490.     if (iszero(z2))
  491.         error("Division by zero");
  492.     if (iszero(z1)) {
  493.         *res = _zero_;
  494.         *rem = _zero_;
  495.         return;
  496.     }
  497.     if (isone(z2)) {
  498.         zcopy(z1, res);
  499.         *rem = _zero_;
  500.         return;
  501.     }
  502.     i = 32768;
  503.     j = 0;
  504.     k = (long) z2.v[z2.len - 1];
  505.     while (! (k & i)) {
  506.         j ++;
  507.         i >>= 1;
  508.     }
  509.     ztmp1.v = alloc(z1.len + 1);
  510.     ztmp1.len = z1.len + 1;
  511.     copyval(z1, ztmp1);
  512.     ztmp1.v[z1.len] = 0;
  513.     ztmp1.sign = 0;
  514.     ztmp2.v = alloc(z2.len);
  515.     ztmp2.len = z2.len;
  516.     ztmp2.sign = 0;
  517.     copyval(z2, ztmp2);
  518.     if (zrel(ztmp1, ztmp2) < 0) {
  519.         rem->v = ztmp1.v;
  520.         rem->sign = z1.sign;
  521.         rem->len = z1.len;
  522.         freeh(ztmp2.v);
  523.         *res = _zero_;
  524.         return;
  525.     }
  526.     quo.len = z1.len - z2.len + 1;
  527.     quo.v = alloc(quo.len);
  528.     quo.sign = z1.sign != z2.sign;
  529.     clearval(quo);
  530.  
  531.     ztmp3.v = zalloctemp(z2.len + 1);
  532.  
  533.     /*
  534.      * Normalize z1 and z2
  535.      */
  536.     shiftl(ztmp1, j);
  537.     shiftl(ztmp2, j);
  538.  
  539.     k = ztmp1.len - ztmp2.len;
  540.     q = quo.v + quo.len;
  541.     y = ztmp1.len - 1;
  542.     h2 = ztmp2.v [ztmp2.len - 1];
  543.     v2 = 0;
  544.     if (ztmp2.len >= 2)
  545.         v2 = ztmp2.v [ztmp2.len - 2];
  546.     for (;k--; --y) {
  547.         pp = ztmp1.v + y - 1;
  548.         pair.silow = pp[0];
  549.         pair.sihigh = pp[1];
  550.         if (ztmp1.v[y] == h2)
  551.             x = BASE1;
  552.         else
  553.             x = pair.ivalue / h2;
  554.         if (x) {
  555.             while (pair.ivalue - x * h2 < BASE &&
  556.                 v2 * x > (pair.ivalue - x * h2) * BASE + ztmp1.v [y-2]) {
  557.                     --x;
  558.             }
  559.             dmul(ztmp2, x, &ztmp3);
  560. #ifdef divblab
  561.             printf(" x: %ld\n", x);
  562.             printf("ztmp1: ");
  563.             printz(ztmp1);
  564.             printf("ztmp2: ");
  565.             printz(ztmp2);
  566.             printf("ztmp3: ");
  567.             printz(ztmp3);
  568. #endif
  569.             if (dsub(ztmp1, ztmp3, y, ztmp2.len)) {
  570.                 --x;
  571.                 /*
  572.                 printf("adding back\n");
  573.                 */
  574.                 dadd(ztmp1, ztmp2, y, ztmp2.len);
  575.             }
  576.         }
  577.         trim(&ztmp1);
  578.         *--q = (HALF)x;
  579.     }
  580.     shiftr(ztmp1, j);
  581.     *rem = ztmp1;
  582.     trim(rem);
  583.     freeh(ztmp2.v);
  584.     trim(&quo);
  585.     *res = quo;
  586. }
  587.  
  588.  
  589. /*
  590.  * Return the quotient and remainder of an integer divided by a small
  591.  * number.  A nonzero remainder is only meaningful when both numbers
  592.  * are positive.
  593.  */
  594. long
  595. zdivi(z, n, res)
  596.     ZVALUE z, *res;
  597.     long n;
  598. {
  599.     register HALF *h1, *sd;
  600.     FULL val;
  601.     HALF divval[2];
  602.     ZVALUE div;
  603.     ZVALUE dest;
  604.     long len;
  605.  
  606.     if (n == 0)
  607.         error("Division by zero");
  608.     if (iszero(z)) {
  609.         *res = _zero_;
  610.         return 0;
  611.     }
  612.     if (n < 0) {
  613.         n = -n;
  614.         z.sign = !z.sign;
  615.     }
  616.     if (n == 1) {
  617.         zcopy(z, res);
  618.         return 0;
  619.     }
  620.     /*
  621.      * If the division is by a large number, then call the normal
  622.      * divide routine.
  623.      */
  624.     if (n & ~BASE1) {
  625.         div.sign = 0;
  626.         div.len = 2;
  627.         div.v = divval;
  628.         divval[0] = (HALF) n;
  629.         divval[1] = ((FULL) n) >> BASEB;
  630.         zdiv(z, div, res, &dest);
  631.         n = (istiny(dest) ? z1tol(dest) : z2tol(dest));
  632.         freeh(dest.v);
  633.         return n;
  634.     }
  635.     /*
  636.      * Division is by a small number, so we can be quick about it.
  637.      */
  638.     len = z.len;
  639.     dest.sign = z.sign;
  640.     dest.len = len;
  641.     dest.v = alloc(len);
  642.     h1 = z.v + len - 1;
  643.     sd = dest.v + len - 1;
  644.     val = 0;
  645.     while (len--) {
  646.         val = ((val << BASEB) + ((FULL) *h1--));
  647.         *sd-- = val / n;
  648.         val %= n;
  649.     }
  650.     quicktrim(dest);
  651.     *res = dest;
  652.     return val;
  653. }
  654.  
  655.  
  656. /*
  657.  * Return the quotient of two numbers.
  658.  * This works the same as zdiv, except that the remainer is not returned.
  659.  */
  660. void
  661. zquo(z1, z2, res)
  662.     ZVALUE z1, z2, *res;
  663. {
  664.     long i, j, k;
  665.     register HALF *q, *pp;
  666.     SIUNION pair;            /* pair of halfword values */
  667.     HALF h2, v2;
  668.     long y;
  669.     FULL x;
  670.     ZVALUE ztmp1, ztmp2, ztmp3, quo;
  671.  
  672.     if (iszero(z2))
  673.         error("Division by zero");
  674.     if (iszero(z1)) {
  675.         *res = _zero_;
  676.         return;
  677.     }
  678.     if (isone(z2)) {
  679.         zcopy(z1, res);
  680.         return;
  681.     }
  682.     i = 32768;
  683.     j = 0;
  684.     k = (long) z2.v[z2.len - 1];
  685.     while (! (k & i)) {
  686.         j ++;
  687.         i >>= 1;
  688.     }
  689.     ztmp1.v = alloc(z1.len + 1);
  690.     ztmp1.len = z1.len + 1;
  691.     copyval(z1, ztmp1);
  692.     ztmp1.v[z1.len] = 0;
  693.     ztmp1.sign = 0;
  694.     ztmp2.v = alloc(z2.len);
  695.     ztmp2.len = z2.len;
  696.     ztmp2.sign = 0;
  697.     copyval(z2, ztmp2);
  698.     if (zrel(ztmp1, ztmp2) < 0) {
  699.         freeh(ztmp1.v);
  700.         freeh(ztmp2.v);
  701.         *res = _zero_;
  702.         return;
  703.     }
  704.     quo.len = z1.len - z2.len + 1;
  705.     quo.v = alloc(quo.len);
  706.     quo.sign = z1.sign != z2.sign;
  707.     clearval(quo);
  708.  
  709.     ztmp3.v = zalloctemp(z2.len + 1);
  710.  
  711.     /*
  712.      * Normalize z1 and z2
  713.      */
  714.     shiftl(ztmp1, j);
  715.     shiftl(ztmp2, j);
  716.  
  717.     k = ztmp1.len - ztmp2.len;
  718.     q = quo.v + quo.len;
  719.     y = ztmp1.len - 1;
  720.     h2 = ztmp2.v [ztmp2.len - 1];
  721.     v2 = 0;
  722.     if (ztmp2.len >= 2)
  723.         v2 = ztmp2.v [ztmp2.len - 2];
  724.     for (;k--; --y) {
  725.         pp = ztmp1.v + y - 1;
  726.         pair.silow = pp[0];
  727.         pair.sihigh = pp[1];
  728.         if (ztmp1.v[y] == h2)
  729.             x = BASE1;
  730.         else
  731.             x = pair.ivalue / h2;
  732.         if (x) {
  733.             while (pair.ivalue - x * h2 < BASE &&
  734.                 v2 * x > (pair.ivalue - x * h2) * BASE + ztmp1.v [y-2]) {
  735.                     --x;
  736.             }
  737.             dmul(ztmp2, x, &ztmp3);
  738.             if (dsub(ztmp1, ztmp3, y, ztmp2.len)) {
  739.                 --x;
  740.                 dadd(ztmp1, ztmp2, y, ztmp2.len);
  741.             }
  742.         }
  743.         trim(&ztmp1);
  744.         *--q = (HALF)x;
  745.     }
  746.     freeh(ztmp1.v);
  747.     freeh(ztmp2.v);
  748.     trim(&quo);
  749.     *res = quo;
  750. }
  751.  
  752.  
  753. /*
  754.  * Compute the remainder after dividing one number by another.
  755.  * This is only defined for positive z2 values.
  756.  * The result is normalized to lie in the range 0 to z2-1.
  757.  */
  758. void
  759. zmod(z1, z2, rem)
  760.     ZVALUE z1, z2, *rem;
  761. {
  762.     long i, j, k, neg;
  763.     register HALF *pp;
  764.     SIUNION pair;            /* pair of halfword values */
  765.     HALF h2, v2;
  766.     long y;
  767.     FULL x;
  768.     ZVALUE ztmp1, ztmp2, ztmp3;
  769.  
  770.     if (iszero(z2))
  771.         error("Division by zero");
  772.     if (isneg(z2))
  773.         error("Non-positive modulus");
  774.     if (iszero(z1) || isunit(z2)) {
  775.         *rem = _zero_;
  776.         return;
  777.     }
  778.     if (istwo(z2)) {
  779.         if (isodd(z1))
  780.             *rem = _one_;
  781.         else
  782.             *rem = _zero_;
  783.         return;
  784.     }
  785.     neg = z1.sign;
  786.     z1.sign = 0;
  787.  
  788.     /*
  789.      * Do a quick check to see if the absolute value of the number
  790.      * is less than the modulus.  If so, then the result is just a
  791.      * subtract or a copy.
  792.      */
  793.     h2 = z1.v[z1.len - 1];
  794.     v2 = z2.v[z2.len - 1];
  795.     if ((z1.len < z2.len) || ((z1.len == z2.len) && (h2 < v2))) {
  796.         if (neg)
  797.             zsub(z2, z1, rem);
  798.         else
  799.             zcopy(z1, rem);
  800.         return;
  801.     }
  802.  
  803.     /*
  804.      * Do another quick check to see if the number is positive and
  805.      * between the size of the modulus and twice the modulus.
  806.      * If so, then the answer is just another subtract.
  807.      */
  808.     if (!neg && (z1.len == z2.len) && (h2 > v2) &&
  809.         (((FULL) h2) < 2 * ((FULL) v2)))
  810.     {
  811.         zsub(z1, z2, rem);
  812.         return;
  813.     }
  814.  
  815.     /*
  816.      * If the modulus is an exact power of two, then the result
  817.      * can be obtained by ignoring the high bits of the number.
  818.      * This truncation assumes that the number of words for the
  819.      * number is at least as large as the number of words in the
  820.      * modulus, which is true at this point.
  821.      */
  822.     if (((v2 & -v2) == v2) && zisonebit(z2)) {    /* ASSUMES 2'S COMP */
  823.         i = zhighbit(z2);
  824.         z1.len = (i + BASEB - 1) / BASEB;
  825.         zcopy(z1, &ztmp1);
  826.         i %= BASEB;
  827.         if (i)
  828.             ztmp1.v[ztmp1.len - 1] &= ((((HALF) 1) << i) - 1);
  829.         ztmp2.len = 0;
  830.         goto gotanswer;
  831.     }
  832.  
  833.     /*
  834.      * If the modulus is one less than an exact power of two, then
  835.      * the result can be simplified similarly to "casting out 9's".
  836.      * Only do this simplification for large enough modulos.
  837.      */
  838.     if ((z2.len > 1) && (z2.v[0] == BASE1) && zisallbits(z2)) {
  839.         i = -(zhighbit(z2) + 1);
  840.         zcopy(z1, &ztmp1);
  841.         z1 = ztmp1;
  842.         while ((k = zrel(z1, z2)) > 0) {
  843.             ztmp1 = _zero_;
  844.             while (!iszero(z1)) {
  845.                 zand(z1, z2, &ztmp2);
  846.                 zadd(ztmp2, ztmp1, &ztmp3);
  847.                 freeh(ztmp1.v);
  848.                 freeh(ztmp2.v);
  849.                 ztmp1 = ztmp3;
  850.                 zshift(z1, i, &ztmp2);
  851.                 freeh(z1.v);
  852.                 z1 = ztmp2;
  853.             }
  854.             freeh(z1.v);
  855.             z1 = ztmp1;
  856.         }
  857.         if (k == 0) {
  858.             freeh(ztmp1.v);
  859.             *rem = _zero_;
  860.             return;
  861.         }
  862.         ztmp2.len = 0;
  863.         goto gotanswer;
  864.     }
  865.  
  866.     /*
  867.      * Must actually do the divide.
  868.      */
  869.     i = 32768;
  870.     j = 0;
  871.     k = (long) z2.v[z2.len - 1];
  872.     while (! (k & i)) {
  873.         j ++;
  874.         i >>= 1;
  875.     }
  876.     ztmp1.v = alloc(z1.len + 1);
  877.     ztmp1.len = z1.len + 1;
  878.     copyval(z1, ztmp1);
  879.     ztmp1.v[z1.len] = 0;
  880.     ztmp1.sign = 0;
  881.     ztmp2.v = alloc(z2.len);
  882.     ztmp2.len = z2.len;
  883.     ztmp2.sign = 0;
  884.     copyval(z2, ztmp2);
  885.     if (zrel(ztmp1, ztmp2) < 0)
  886.         goto gotanswer;
  887.  
  888.     ztmp3.v = zalloctemp(z2.len + 1);
  889.  
  890.     /*
  891.      * Normalize z1 and z2
  892.      */
  893.     shiftl(ztmp1, j);
  894.     shiftl(ztmp2, j);
  895.  
  896.     k = ztmp1.len - ztmp2.len;
  897.     y = ztmp1.len - 1;
  898.     h2 = ztmp2.v [ztmp2.len - 1];
  899.     v2 = 0;
  900.     if (ztmp2.len >= 2)
  901.         v2 = ztmp2.v [ztmp2.len - 2];
  902.     for (;k--; --y) {
  903.         pp = ztmp1.v + y - 1;
  904.         pair.silow = pp[0];
  905.         pair.sihigh = pp[1];
  906.         if (ztmp1.v[y] == h2)
  907.             x = BASE1;
  908.         else
  909.             x = pair.ivalue / h2;
  910.         if (x) {
  911.             while (pair.ivalue - x * h2 < BASE &&
  912.                 v2 * x > (pair.ivalue - x * h2) * BASE + ztmp1.v [y-2]) {
  913.                     --x;
  914.             }
  915.             dmul(ztmp2, x, &ztmp3);
  916.             if (dsub(ztmp1, ztmp3, y, ztmp2.len))
  917.                 dadd(ztmp1, ztmp2, y, ztmp2.len);
  918.         }
  919.         trim(&ztmp1);
  920.     }
  921.     shiftr(ztmp1, j);
  922.  
  923. gotanswer:
  924.     trim(&ztmp1);
  925.     if (ztmp2.len)
  926.         freeh(ztmp2.v);
  927.     if (neg && !iszero(ztmp1)) {
  928.         zsub(z2, ztmp1, rem);
  929.         freeh(ztmp1.v);
  930.     } else
  931.         *rem = ztmp1;
  932. }
  933.  
  934.  
  935. /*
  936.  * Calculate the mod of an integer by a small number.
  937.  * This is only defined for positive moduli.
  938.  */
  939. long
  940. zmodi(z, n)
  941.     ZVALUE z;
  942.     long n;
  943. {
  944.     register HALF *h1;
  945.     FULL val;
  946.     HALF divval[2];
  947.     ZVALUE div;
  948.     ZVALUE temp;
  949.     long len;
  950.  
  951.     if (n == 0)
  952.         error("Division by zero");
  953.     if (n < 0)
  954.         error("Non-positive modulus");
  955.     if (iszero(z) || (n == 1))
  956.         return 0;
  957.     if (isone(z))
  958.         return 1;
  959.     /*
  960.      * If the modulus is by a large number, then call the normal
  961.      * modulo routine.
  962.      */
  963.     if (n & ~BASE1) {
  964.         div.sign = 0;
  965.         div.len = 2;
  966.         div.v = divval;
  967.         divval[0] = (HALF) n;
  968.         divval[1] = ((FULL) n) >> BASEB;
  969.         zmod(z, div, &temp);
  970.         n = (istiny(temp) ? z1tol(temp) : z2tol(temp));
  971.         freeh(temp.v);
  972.         return n;
  973.     }
  974.     /*
  975.      * The modulus is by a small number, so we can do this quickly.
  976.      */
  977.     len = z.len;
  978.     h1 = z.v + len - 1;
  979.     val = 0;
  980.     while (len--)
  981.         val = ((val << BASEB) + ((FULL) *h1--)) % n;
  982.     if (z.sign)
  983.         val = n - val;
  984.     return val;
  985. }
  986.  
  987.  
  988. /*
  989.  * Return whether or not one number exactly divides another one.
  990.  * Returns TRUE if division occurs with no remainder.
  991.  * z1 is the number to be divided by z2.
  992.  */
  993. BOOL
  994. zdivides(z1, z2)
  995.     ZVALUE z1, z2;        /* numbers to test division into and by */
  996. {
  997.     ZVALUE temp;
  998.     long cv;
  999.  
  1000.     z1.sign = 0;
  1001.     z2.sign = 0;
  1002.     /*
  1003.      * Take care of obvious cases first
  1004.      */
  1005.     if (isleone(z2)) {    /* division by zero or one */
  1006.         if (*z2.v == 0)
  1007.             error("Division by zero");
  1008.         return TRUE;
  1009.     }
  1010.     if (iszero(z1))    /* everything divides zero */
  1011.         return TRUE;
  1012.     if (z1.len < z2.len)    /* quick size comparison */
  1013.         return FALSE;
  1014.     if ((z1.len == z2.len) && (z1.v[z1.len-1] < z2.v[z2.len-1]))    /* more */
  1015.         return FALSE;
  1016.     if (isodd(z1) && iseven(z2))    /* can't divide odd by even */
  1017.         return FALSE;
  1018.     if (zlowbit(z1) < zlowbit(z2))    /* can't have smaller power of two */
  1019.         return FALSE;
  1020.     cv = zrel(z1, z2);    /* can't divide smaller number */
  1021.     if (cv <= 0)
  1022.         return (cv == 0);
  1023.     /*
  1024.      * Now do the real work.  Divisor divides dividend if the gcd of the
  1025.      * two numbers equals the divisor.
  1026.      */
  1027.     zgcd(z1, z2, &temp);
  1028.     cv = zcmp(z2, temp);
  1029.     freeh(temp.v);
  1030.     return (cv == 0);
  1031. }
  1032.  
  1033.  
  1034. /*
  1035.  * Compute the logical OR of two numbers
  1036.  */
  1037. void
  1038. zor(z1, z2, res)
  1039.     ZVALUE z1, z2, *res;
  1040. {
  1041.     register HALF *sp, *dp;
  1042.     long len;
  1043.     ZVALUE bz, lz, dest;
  1044.  
  1045.     if (z1.len >= z2.len) {
  1046.         bz = z1;
  1047.         lz = z2;
  1048.     } else {
  1049.         bz = z2;
  1050.         lz = z1;
  1051.     }
  1052.     dest.len = bz.len;
  1053.     dest.v = alloc(dest.len);
  1054.     dest.sign = 0;
  1055.     copyval(bz, dest);
  1056.     len = lz.len;
  1057.     sp = lz.v;
  1058.     dp = dest.v;
  1059.     while (len--)
  1060.         *dp++ |= *sp++;
  1061.     *res = dest;
  1062. }
  1063.  
  1064.  
  1065. /*
  1066.  * Compute the logical AND of two numbers.
  1067.  */
  1068. void
  1069. zand(z1, z2, res)
  1070.     ZVALUE z1, z2, *res;
  1071. {
  1072.     register HALF *h1, *h2, *hd;
  1073.     register long len;
  1074.     ZVALUE dest;
  1075.  
  1076.     len = ((z1.len <= z2.len) ? z1.len : z2.len);
  1077.     h1 = &z1.v[len-1];
  1078.     h2 = &z2.v[len-1];
  1079.     while ((len > 1) && ((*h1 & *h2) == 0)) {
  1080.         h1--;
  1081.         h2--;
  1082.         len--;
  1083.     }
  1084.     dest.len = len;
  1085.     dest.v = alloc(len);
  1086.     dest.sign = 0;
  1087.     h1 = z1.v;
  1088.     h2 = z2.v;
  1089.     hd = dest.v;
  1090.     while (len--)
  1091.         *hd++ = (*h1++ & *h2++);
  1092.     *res = dest;
  1093. }
  1094.  
  1095.  
  1096. /*
  1097.  * Compute the logical XOR of two numbers.
  1098.  */
  1099. void
  1100. zxor(z1, z2, res)
  1101.     ZVALUE z1, z2, *res;
  1102. {
  1103.     register HALF *sp, *dp;
  1104.     long len;
  1105.     ZVALUE bz, lz, dest;
  1106.  
  1107.     if (z1.len == z2.len) {
  1108.         for (len = z1.len; ((len > 1) && (z1.v[len-1] == z2.v[len-1])); len--) ;
  1109.         z1.len = len;
  1110.         z2.len = len;
  1111.     }
  1112.     if (z1.len >= z2.len) {
  1113.         bz = z1;
  1114.         lz = z2;
  1115.     } else {
  1116.         bz = z2;
  1117.         lz = z1;
  1118.     }
  1119.     dest.len = bz.len;
  1120.     dest.v = alloc(dest.len);
  1121.     dest.sign = 0;
  1122.     copyval(bz, dest);
  1123.     len = lz.len;
  1124.     sp = lz.v;
  1125.     dp = dest.v;
  1126.     while (len--)
  1127.         *dp++ ^= *sp++;
  1128.     *res = dest;
  1129. }
  1130.  
  1131.  
  1132. /*
  1133.  * Shift a number left (or right) by the specified number of bits.
  1134.  * Positive shift means to the left.  When shifting right, rightmost
  1135.  * bits are lost.  The sign of the number is preserved.
  1136.  */
  1137. void
  1138. zshift(z, n, res)
  1139.     ZVALUE z, *res;
  1140.     long n;
  1141. {
  1142.     ZVALUE ans;
  1143.     long hc;        /* number of halfwords shift is by */
  1144.  
  1145.     if (iszero(z)) {
  1146.         *res = _zero_;
  1147.         return;
  1148.     }
  1149.     if (n == 0) {
  1150.         zcopy(z, res);
  1151.         return;
  1152.     }
  1153.     /*
  1154.      * If shift value is negative, then shift right.
  1155.      * Check for large shifts, and handle word-sized shifts quickly.
  1156.      */
  1157.     if (n < 0) {
  1158.         n = -n;
  1159.         if ((n < 0) || (n >= (z.len * BASEB))) {
  1160.             *res = _zero_;
  1161.             return;
  1162.         }
  1163.         hc = n / BASEB;
  1164.         n %= BASEB;
  1165.         z.v += hc;
  1166.         z.len -= hc;
  1167.         ans.len = z.len;
  1168.         ans.v = alloc(ans.len);
  1169.         ans.sign = z.sign;
  1170.         copyval(z, ans);
  1171.         if (n > 0) {
  1172.             shiftr(ans, n);
  1173.             trim(&ans);
  1174.         }
  1175.         if (iszero(ans)) {
  1176.             freeh(ans.v);
  1177.             ans = _zero_;
  1178.         }
  1179.         *res = ans;
  1180.         return;
  1181.     }
  1182.     /*
  1183.      * Shift value is positive, so shift leftwards.
  1184.      * Check specially for a shift of the value 1, since this is common.
  1185.      * Also handle word-sized shifts quickly.
  1186.      */
  1187.     if (isunit(z)) {
  1188.         zbitvalue(n, res);
  1189.         res->sign = z.sign;
  1190.         return;
  1191.     }
  1192.     hc = n / BASEB;
  1193.     n %= BASEB;
  1194.     ans.len = z.len + hc + 1;
  1195.     ans.v = alloc(ans.len);
  1196.     ans.sign = z.sign;
  1197.     if (hc > 0)
  1198.         memset((char *) ans.v, 0, hc * sizeof(HALF));
  1199.     memcpy((char *) (ans.v + hc), 
  1200.         (char *) z.v, z.len * sizeof(HALF));
  1201.     ans.v[ans.len - 1] = 0;
  1202.     if (n > 0) {
  1203.         ans.v += hc;
  1204.         ans.len -= hc;
  1205.         shiftl(ans, n);
  1206.         ans.v -= hc;
  1207.         ans.len += hc;
  1208.     }
  1209.     trim(&ans);
  1210.     *res = ans;
  1211. }
  1212.  
  1213.  
  1214. /*
  1215.  * Return the position of the lowest bit which is set in the binary
  1216.  * representation of a number (counting from zero).  This is the highest
  1217.  * power of two which evenly divides the number.
  1218.  */
  1219. long
  1220. zlowbit(z)
  1221.     ZVALUE z;
  1222. {
  1223.     register HALF *zp;
  1224.     long n;
  1225.     HALF dataval;
  1226.     HALF *bitval;
  1227.  
  1228.     n = 0;
  1229.     for (zp = z.v; *zp == 0; zp++)
  1230.         if (++n >= z.len)
  1231.             return 0;
  1232.     dataval = *zp;
  1233.     bitval = bitmask;
  1234.     while ((*(bitval++) & dataval) == 0) {
  1235.     }
  1236.     return (n*BASEB)+(bitval-bitmask-1);
  1237. }
  1238.  
  1239.  
  1240. /*
  1241.  * Return the position of the highest bit which is set in the binary
  1242.  * representation of a number (counting from zero).  This is the highest power
  1243.  * of two which is less than or equal to the number (which is assumed nonzero).
  1244.  */
  1245. long
  1246. zhighbit(z)
  1247.     ZVALUE z;
  1248. {
  1249.     HALF dataval;
  1250.     HALF *bitval;
  1251.  
  1252.     dataval = z.v[z.len-1];
  1253.     if (dataval) {
  1254.         bitval = bitmask+BASEB;
  1255.         while ((*(--bitval) & dataval) == 0) {
  1256.         }
  1257.     }
  1258.     return (z.len*BASEB)+(bitval-bitmask-BASEB);
  1259. }
  1260.  
  1261.  
  1262. #if 0
  1263. /*
  1264.  * Reverse the bits of a particular range of bits of a number.
  1265.  *
  1266.  * This function returns an integer with bits a thru b swapped.
  1267.  * That is, bit a is swapped with bit b, bit a+1 is swapped with b-1,
  1268.  * and so on.
  1269.  *
  1270.  * As a special case, if the ending bit position is < 0, is it taken to 
  1271.  * mean the highest bit set.  Thus zbitrev(0, -1, z, &res) will 
  1272.  * perform a complete bit reverse of the number 'z'.
  1273.  *
  1274.  * As a special case, if the starting bit position is < 0, is it taken to 
  1275.  * mean the lowest bit set.  Thus zbitrev(-1, -1, z, &res) is the
  1276.  * same as zbitrev(lowbit(z), highbit(z), z, &res).
  1277.  *
  1278.  * Note that the low order bit number is taken to be 0.  Also, bitrev
  1279.  * ignores the sign of the number.
  1280.  *
  1281.  * Bits beyond the highest bit are taken to be zero.  Thus the calling
  1282.  * bitrev(0, 100, _one_, &res) will result in a value of 2^100.
  1283.  */
  1284. void
  1285. zbitrev(low, high, z, res)
  1286.     long low;    /* lowest bit to reverse, <0 => lowbit(z) */
  1287.     long high;    /* highest bit to reverse, <0 => highbit(z) */
  1288.     ZVALUE z;    /* value to bit reverse */
  1289.     ZVALUE *res;    /* resulting bit reverse number */
  1290. {
  1291. }
  1292. #endif
  1293.  
  1294.  
  1295. /*
  1296.  * Return whether or not the specifed bit number is set in a number.
  1297.  * Rightmost bit of a number is bit 0.
  1298.  */
  1299. BOOL
  1300. zisset(z, n)
  1301.     ZVALUE z;
  1302.     long n;
  1303. {
  1304.     if ((n < 0) || ((n / BASEB) >= z.len))
  1305.         return FALSE;
  1306.     return ((z.v[n / BASEB] & (((HALF) 1) << (n % BASEB))) != 0);
  1307. }
  1308.  
  1309.  
  1310. /*
  1311.  * Check whether or not a number has exactly one bit set, and
  1312.  * thus is an exact power of two.  Returns TRUE if so.
  1313.  */
  1314. BOOL
  1315. zisonebit(z)
  1316.     ZVALUE z;
  1317. {
  1318.     register HALF *hp;
  1319.     register LEN len;
  1320.  
  1321.     if (iszero(z) || isneg(z))
  1322.         return FALSE;
  1323.     hp = z.v;
  1324.     len = z.len;
  1325.     while (len > 4) {
  1326.         len -= 4;
  1327.         if (*hp++ || *hp++ || *hp++ || *hp++)
  1328.             return FALSE;
  1329.     }
  1330.     while (--len > 0) {
  1331.         if (*hp++)
  1332.             return FALSE;
  1333.     }
  1334.     return ((*hp & -*hp) == *hp);        /* NEEDS 2'S COMPLEMENT */
  1335. }
  1336.  
  1337.  
  1338. /*
  1339.  * Check whether or not a number has all of its bits set below some
  1340.  * bit position, and thus is one less than an exact power of two.
  1341.  * Returns TRUE if so.
  1342.  */
  1343. BOOL
  1344. zisallbits(z)
  1345.     ZVALUE z;
  1346. {
  1347.     register HALF *hp;
  1348.     register LEN len;
  1349.     HALF digit;
  1350.  
  1351.     if (iszero(z) || isneg(z))
  1352.         return FALSE;
  1353.     hp = z.v;
  1354.     len = z.len;
  1355.     while (len > 4) {
  1356.         len -= 4;
  1357.         if ((*hp++ != BASE1) || (*hp++ != BASE1) ||
  1358.             (*hp++ != BASE1) || (*hp++ != BASE1))
  1359.                 return FALSE;
  1360.     }
  1361.     while (--len > 0) {
  1362.         if (*hp++ != BASE1)
  1363.             return FALSE;
  1364.     }
  1365.     digit = *hp + 1;
  1366.     return ((digit & -digit) == digit);    /* NEEDS 2'S COMPLEMENT */
  1367. }
  1368.  
  1369.  
  1370. /*
  1371.  * Return the number whose binary representation contains only one bit which
  1372.  * is in the specified position (counting from zero).  This is equivilant
  1373.  * to raising two to the given power.
  1374.  */
  1375. void
  1376. zbitvalue(n, res)
  1377.     long n;
  1378.     ZVALUE *res;
  1379. {
  1380.     ZVALUE z;
  1381.  
  1382.     if (n < 0) n = 0;
  1383.     z.sign = 0;
  1384.     z.len = (n / BASEB) + 1;
  1385.     z.v = alloc(z.len);
  1386.     clearval(z);
  1387.     z.v[z.len-1] = (((HALF) 1) << (n % BASEB));
  1388.     *res = z;
  1389. }
  1390.  
  1391.  
  1392. /*
  1393.  * Compare a number against zero.
  1394.  * Returns the sgn function of the number (-1, 0, or 1).
  1395.  */
  1396. FLAG
  1397. ztest(z)
  1398.     ZVALUE z;
  1399. {
  1400.     register int sign;
  1401.     register HALF *h;
  1402.     register long len;
  1403.  
  1404.     sign = 1;
  1405.     if (z.sign)
  1406.         sign = -sign;
  1407.     h = z.v;
  1408.     len = z.len;
  1409.     while (len--)
  1410.         if (*h++)
  1411.             return sign;
  1412.     return 0;
  1413. }
  1414.  
  1415.  
  1416. /*
  1417.  * Compare two numbers to see which is larger.
  1418.  * Returns -1 if first number is smaller, 0 if they are equal, and 1 if
  1419.  * first number is larger.  This is the same result as ztest(z2-z1).
  1420.  */
  1421. FLAG
  1422. zrel(z1, z2)
  1423.     ZVALUE z1, z2;
  1424. {
  1425.     register HALF *h1, *h2;
  1426.     register long len1, len2;
  1427.     int sign;
  1428.  
  1429.     sign = 1;
  1430.     if (z1.sign < z2.sign)
  1431.         return 1;
  1432.     if (z2.sign < z1.sign)
  1433.         return -1;
  1434.     if (z2.sign)
  1435.         sign = -1;
  1436.     len1 = z1.len;
  1437.     len2 = z2.len;
  1438.     h1 = z1.v + z1.len - 1;
  1439.     h2 = z2.v + z2.len - 1;
  1440.     while (len1 > len2) {
  1441.         if (*h1--)
  1442.             return sign;
  1443.         len1--;
  1444.     }
  1445.     while (len2 > len1) {
  1446.         if (*h2--)
  1447.             return -sign;
  1448.         len2--;
  1449.     }
  1450.     while (len1--) {
  1451.         if (*h1-- != *h2--)
  1452.             break;
  1453.     }
  1454.     if ((len1 = *++h1) > (len2 = *++h2))
  1455.         return sign;
  1456.     if (len1 < len2)
  1457.         return -sign;
  1458.     return 0;
  1459. }
  1460.  
  1461.  
  1462. /*
  1463.  * Compare two numbers to see if they are equal or not.
  1464.  * Returns TRUE if they differ.
  1465.  */
  1466. BOOL
  1467. zcmp(z1, z2)
  1468.     ZVALUE z1, z2;
  1469. {
  1470.     register HALF *h1, *h2;
  1471.     register long len;
  1472.  
  1473.     if ((z1.sign != z2.sign) || (z1.len != z2.len) || (*z1.v != *z2.v))
  1474.         return TRUE;
  1475.     len = z1.len;
  1476.     h1 = z1.v;
  1477.     h2 = z2.v;
  1478.     while (len-- > 0) {
  1479.         if (*h1++ != *h2++)
  1480.             return TRUE;
  1481.     }
  1482.     return FALSE;
  1483. }
  1484.  
  1485.  
  1486. /*
  1487.  * Internal utility subroutines
  1488.  */
  1489. static void
  1490. dadd(z1, z2, y, n)
  1491.     ZVALUE z1, z2;
  1492.     long y, n;
  1493. {
  1494.     HALF *s1, *s2;
  1495.     short carry;
  1496.     long sum;
  1497.  
  1498.     s1 = z1.v + y - n;
  1499.     s2 = z2.v;
  1500.     carry = 0;
  1501.     while (n--) {
  1502.         sum = (long)*s1 + (long)*s2 + carry;
  1503.         carry = 0;
  1504.         if (sum >= BASE) {
  1505.             sum -= BASE;
  1506.             carry = 1;
  1507.         }
  1508.         *s1 = (HALF)sum;
  1509.         ++s1;
  1510.         ++s2;
  1511.     }
  1512.     sum = (long)*s1 + carry;
  1513.     *s1 = (HALF)sum;
  1514. }
  1515.  
  1516.  
  1517. /*
  1518.  * Do subtract for divide, returning TRUE if subtraction went negative.
  1519.  */
  1520. static BOOL
  1521. dsub(z1, z2, y, n)
  1522.     ZVALUE z1, z2;
  1523.     long y, n;
  1524. {
  1525.     HALF *s1, *s2, *s3;
  1526.     FULL i1;
  1527.     BOOL neg;
  1528.  
  1529.     neg = FALSE;
  1530.     s1 = z1.v + y - n;
  1531.     s2 = z2.v;
  1532.     if (++n > z2.len)
  1533.         n = z2.len;
  1534.     while (n--) {
  1535.         i1 = (FULL) *s1;
  1536.         if (i1 < (FULL) *s2) {
  1537.             s3 = s1 + 1;
  1538.             while (s3 < z1.v + z1.len && !(*s3)) {
  1539.                 *s3 = BASE1;
  1540.                 ++s3;
  1541.             }
  1542.             if (s3 >= z1.v + z1.len)
  1543.                 neg = TRUE;
  1544.             else
  1545.                 --(*s3);
  1546.             i1 += BASE;
  1547.         }
  1548.         *s1 = i1 - (FULL) *s2;
  1549.         ++s1;
  1550.         ++s2;
  1551.     }
  1552.     return neg;
  1553. }
  1554.  
  1555.  
  1556. /*
  1557.  * Multiply a number by a single 'digit'.
  1558.  * This is meant to be used only by the divide routine, and so the
  1559.  * destination area must already be allocated and be large enough.
  1560.  */
  1561. static void
  1562. dmul(z, mul, dest)
  1563.     ZVALUE z;
  1564.     FULL mul;
  1565.     ZVALUE *dest;
  1566. {
  1567.     register HALF *zp, *dp;
  1568.     SIUNION sival;
  1569.     FULL carry;
  1570.     long len;
  1571.  
  1572.     dp = dest->v;
  1573.     dest->sign = 0;
  1574.     if (mul == 0) {
  1575.         dest->len = 1;
  1576.         *dp = 0;
  1577.         return;
  1578.     }
  1579.     len = z.len;
  1580.     zp = z.v + len - 1;
  1581.     while ((*zp == 0) && (len > 1)) {
  1582.         len--;
  1583.         zp--;
  1584.     }
  1585.     dest->len = len;
  1586.     zp = z.v;
  1587.     carry = 0;
  1588.     while (len >= 4) {
  1589.         len -= 4;
  1590.         sival.ivalue = (mul * ((FULL) *zp++)) + carry;
  1591.         *dp++ = sival.silow;
  1592.         sival.ivalue = (mul * ((FULL) *zp++)) + ((FULL) sival.sihigh);
  1593.         *dp++ = sival.silow;
  1594.         sival.ivalue = (mul * ((FULL) *zp++)) + ((FULL) sival.sihigh);
  1595.         *dp++ = sival.silow;
  1596.         sival.ivalue = (mul * ((FULL) *zp++)) + ((FULL) sival.sihigh);
  1597.         *dp++ = sival.silow;
  1598.         carry = sival.sihigh;
  1599.     }
  1600.     while (--len >= 0) {
  1601.         sival.ivalue = (mul * ((FULL) *zp++)) + carry;
  1602.         *dp++ = sival.silow;
  1603.         carry = sival.sihigh;
  1604.     }
  1605.     if (carry) {
  1606.         *dp = (HALF)carry;
  1607.         dest->len++;
  1608.     }
  1609. }
  1610.  
  1611.  
  1612. /*
  1613.  * Utility to calculate the gcd of two small integers.
  1614.  */
  1615. long
  1616. iigcd(i1, i2)
  1617.     long i1, i2;
  1618. {
  1619.     FULL f1, f2, temp;
  1620.  
  1621.     f1 = (FULL) ((i1 >= 0) ? i1 : -i1);
  1622.     f2 = (FULL) ((i2 >= 0) ? i2 : -i2);
  1623.     while (f1) {
  1624.         temp = f2 % f1;
  1625.         f2 = f1;
  1626.         f1 = temp;
  1627.     }
  1628.     return (long) f2;
  1629. }
  1630.  
  1631.  
  1632. void
  1633. trim(z)
  1634.     ZVALUE *z;
  1635. {
  1636.     register HALF *h;
  1637.     register long len;
  1638.  
  1639.     h = z->v + z->len - 1;
  1640.     len = z->len;
  1641.     while (*h == 0 && len > 1) {
  1642.         --h;
  1643.         --len;
  1644.     }
  1645.     z->len = len;
  1646. }
  1647.  
  1648.  
  1649. /*
  1650.  * Utility routine to shift right.
  1651.  */
  1652. void
  1653. shiftr(z, n)
  1654.     ZVALUE z;
  1655.     long n;
  1656. {
  1657.     register HALF *h, *lim;
  1658.     FULL mask, maskt;
  1659.     long len;
  1660.  
  1661.     if (n >= BASEB) {
  1662.         len = n / BASEB;
  1663.         h = z.v;
  1664.         lim = z.v + z.len - len;
  1665.         while (h < lim) {
  1666.             h[0] = h[len];
  1667.             ++h;
  1668.         }
  1669.         n -= BASEB * len;
  1670.         lim = z.v + z.len;
  1671.         while (h < lim)
  1672.             *h++ = 0;
  1673.     }
  1674.     if (n) {
  1675.         len = z.len;
  1676.         h = z.v + len - 1;
  1677.         mask = 0;
  1678.         while (len--) {
  1679.             maskt = (((FULL) *h) << (BASEB - n)) & BASE1;
  1680.             *h = (*h >> n) | mask;
  1681.             mask = maskt;
  1682.             --h;
  1683.         }
  1684.     }
  1685. }
  1686.  
  1687.  
  1688. /*
  1689.  * Utility routine to shift left.
  1690.  */
  1691. void
  1692. shiftl(z, n)
  1693.     ZVALUE z;
  1694.     long n;
  1695. {
  1696.     register HALF *h;
  1697.     FULL mask, i;
  1698.     long len;
  1699.  
  1700.     if (n >= BASEB) {
  1701.         len = n / BASEB;
  1702.         h = z.v + z.len - 1;
  1703.         while (!*h)
  1704.             --h;
  1705.         while (h >= z.v) {
  1706.             h[len] = h[0];
  1707.             --h;
  1708.         }
  1709.         n -= BASEB * len;
  1710.         while (len)
  1711.             h[len--] = 0;
  1712.     }
  1713.     if (n > 0) {
  1714.         len = z.len;
  1715.         h = z.v;
  1716.         mask = 0;
  1717.         while (len--) {
  1718.             i = (((FULL) *h) << n) | mask;
  1719.             if (i > BASE1) {
  1720.                 mask = i >> BASEB;
  1721.                 i &= BASE1;
  1722.             } else
  1723.                 mask = 0;
  1724.             *h = (HALF) i;
  1725.             ++h;
  1726.         }
  1727.     }
  1728. }
  1729.  
  1730. /*
  1731.  * initmasks - init the bitmask rotation arrays
  1732.  *
  1733.  * bitmask[i]      (1 << (i-1)),  for  -BASEB*4<=i<=BASEB*4
  1734.  * rotmask[j][k] (1 << ((j+k-1)%BASEB)),  for  -BASEB*2<=j,k<=BASEB*2
  1735.  *
  1736.  * The bmask array contains 8 cycles of rotations of a bit mask.
  1737.  * We point bitmask and the rotmask pointers into the middle to
  1738.  * ensure that we can have a complete two cycle swing forward
  1739.  * and backward.
  1740.  */
  1741. void
  1742. initmasks()
  1743. {
  1744.     int i;
  1745.  
  1746.     /*
  1747.      * setup the bmask array
  1748.      */
  1749.     bmask = alloc((8*BASEB)+1);
  1750.     for (i=0; i < (8*BASEB)+1; ++i) {
  1751.         bmask[i] = 1 << (i%BASEB);
  1752.     }
  1753.  
  1754.     /*
  1755.      * setup the rmask pointers
  1756.      */
  1757.     rmask = (HALF **)malloc(sizeof(HALF *)*((BASEB*4)+2));
  1758.     for (i = 0; i <= (4*BASEB)+1; ++i) {
  1759.         rmask[i] = &bmask[(2*BASEB)+i];
  1760.     }
  1761.  
  1762. #if 0
  1763.     /* 
  1764.      * setup the rotmask array, allow for -2*BASEB thru 2*BASEB indexing
  1765.      */
  1766.     rotmask = &rmask[2*BASEB];
  1767. #endif
  1768.  
  1769.     /*
  1770.      * setup the bitmask array to allow -4*BASEB thru 4*BASEB indexing
  1771.      */
  1772.     bitmask = &bmask[4*BASEB];
  1773.     return;
  1774. }
  1775.  
  1776. /* END CODE */
  1777.